module SF_SP_LINALG
  USE SF_MISC,   only: assert_shape
  USE SF_RANDOM, only: mt_random
  USE SF_LINALG, only: eye,eigh
#ifdef _MPI
  USE SF_MPI
#endif
  implicit none
#ifdef _MPI
  include 'mpif.h'
#endif
  private
  

  interface sp_eigh
     module procedure :: lanczos_arpack_d
     module procedure :: lanczos_arpack_c
#ifdef _MPI
     module procedure :: lanczos_parpack_d
     module procedure :: lanczos_parpack_c
#endif
  end interface sp_eigh


  interface sp_lanc_eigh
     module procedure :: lanczos_eigh_d
     module procedure :: lanczos_eigh_c
#ifdef _MPI
     module procedure :: mpi_lanczos_eigh_d
     module procedure :: mpi_lanczos_eigh_c
#endif
  end interface sp_lanc_eigh


  interface sp_lanc_tridiag
     module procedure :: lanczos_tridiag_d
     module procedure :: lanczos_tridiag_c
#ifdef _MPI
     module procedure :: mpi_lanczos_tridiag_d
     module procedure :: mpi_lanczos_tridiag_c
#endif
  end interface sp_lanc_tridiag


  interface sp_dvdson_eigh
     module procedure :: dvdson_eigh_d
  end interface sp_dvdson_eigh


  complex(8),parameter              :: zero=(0d0,0d0)
  complex(8),parameter              :: one=(1d0,0d0)
  complex(8),parameter              :: xi=(0d0,1d0)
  integer,allocatable               :: seed_random(:)
  integer                           :: nrandom
  logical                           :: verb=.false.
  real(8)                           :: threshold_=1.d-12
  integer                           :: ncheck_=10



  !****************************************************************************************
  !                                      PUBLIC 
  !****************************************************************************************
  public :: sp_eigh
  public :: sp_lanc_eigh
  public :: sp_lanc_tridiag
  public :: sp_dvdson_eigh
  !****************************************************************************************




contains


  !##################################################################
  ! ARPACK METHOD for LOWEST part of the spectrum of a of a SPARSE
  !   MATRIX (defined via H*v)
  ! - DBLE and CMPLX versions included
  ! - SERIAL and PARALLEL-MPI versions included
  !        [COMM    MPI  Communicator for the processor grid.  (INPUT)]
  include "arpack_d.f90"
  include "arpack_c.f90"
#ifdef _MPI
  include "parpack_d.f90"
  include "parpack_c.f90"
#endif


  !##################################################################
  ! LANCZOS METHOD for LOWEST EigenSolution OR tri-diagonalization
  !    of a SPARSE MATRIX (defined via H*v)
  ! - DBLE and CMPLX versions included
  ! - SERIAL and PARALLEL-MPI versions included
  include "lanczos_d.f90"
  include "lanczos_c.f90"
#ifdef _MPI
  include "mpi_lanczos_d.f90"
  include "mpi_lanczos_c.f90"
#endif


  !##################################################################
  ! DAVIDSON METHOD for LOWEST EigenSolution of a SPARSE MATRIX (defined via H*v)
  include "dvdson_serial.f90"



end module SF_SP_LINALG
























! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! !++++++++++++++++++COMPUTATIONAL ROUTINE: TQL2++++++++++++++++++++++++ 
! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! !---------------------------------------------------------------------
! ! PURPOSE computes all eigenvalues/vectors, real symmetric tridiagonal matrix.
! !    This subroutine finds the eigenvalues and eigenvectors of a symmetric
! !    tridiagonal matrix by the QL method.  The eigenvectors of a full
! !    symmetric matrix can also be found if TRED2 has been used to reduce this
! !    full matrix to tridiagonal form.
! !  Parameters:
! !    Input, integer ( kind = 4 ) N, the order of the matrix.
! !
! !    Input/output, real ( kind = 8 ) D(N).  On input, the diagonal elements of
! !    the matrix.  On output, the eigenvalues in ascending order.  If an error
! !    exit is made, the eigenvalues are correct but unordered for indices
! !    1,2,...,IERR-1.
! !
! !    Input/output, real ( kind = 8 ) E(N).  On input, E(2:N) contains the
! !    subdiagonal elements of the input matrix, and E(1) is arbitrary.
! !    On output, E has been destroyed.
! !
! !    Input, real ( kind = 8 ) Z(N,N).  On input, the transformation matrix
! !    produced in the reduction by TRED2, if performed.  If the eigenvectors of
! !    the tridiagonal matrix are desired, Z must contain the identity matrix.
! !    On output, Z contains the orthonormal eigenvectors of the symmetric
! !    tridiagonal (or full) matrix.  If an error exit is made, Z contains
! !    the eigenvectors associated with the stored eigenvalues.
! !
! !    Output, integer ( kind = 4 ) IERR, error flag.
! !    0, normal return,
! !    J, if the J-th eigenvalue has not been determined after
! !    30 iterations.
! !
! !---------------------------------------------------------------------
! subroutine tql2 ( n, d, e, z, ierr )
!   integer :: n
!   real(8) :: c
!   real(8) :: c2
!   real(8) :: c3
!   real(8) :: d(n)
!   real(8) :: dl1
!   real(8) :: e(n)
!   real(8) :: el1
!   real(8) :: f
!   real(8) :: g
!   real(8) :: h
!   integer ( kind = 4 ) i
!   integer ( kind = 4 ) ierr
!   integer ( kind = 4 ) ii
!   integer ( kind = 4 ) j
!   integer ( kind = 4 ) k
!   integer ( kind = 4 ) l
!   integer ( kind = 4 ) l1
!   integer ( kind = 4 ) l2
!   integer ( kind = 4 ) m
!   integer ( kind = 4 ) mml
!   real(8) :: p
!   real(8) :: r
!   real(8) :: s
!   real(8) :: s2
!   real(8) :: tst1
!   real(8) :: tst2
!   real(8) :: z(n,n)
!   ierr = 0
!   if ( n == 1 ) then
!      return
!   end if
!   do i = 2, n
!      e(i-1) = e(i)
!   end do
!   f = 0.0D+00
!   tst1 = 0.0D+00
!   e(n) = 0.0D+00
!   do l = 1, n
!      j = 0
!      h = abs ( d(l) ) + abs ( e(l) )
!      tst1 = max ( tst1, h )
!      !
!      !  Look for a small sub-diagonal element.
!      !
!      do m = l, n
!         tst2 = tst1 + abs ( e(m) )
!         if ( tst2 == tst1 ) then
!            exit
!         end if
!      end do
!      if ( m == l ) then
!         go to 220
!      end if
! 130  continue
!      if ( 30 <= j ) then
!         ierr = l
!         return
!      end if
!      j = j + 1
!      !
!      !  Form shift.
!      !
!      l1 = l + 1
!      l2 = l1 + 1
!      g = d(l)
!      p = ( d(l1) - g ) / ( 2.0D+00 * e(l) )
!      r = pythag ( p, 1.0D+00 )
!      d(l) = e(l) / ( p + sign ( r, p ) )
!      d(l1) = e(l) * ( p + sign ( r, p ) )
!      dl1 = d(l1)
!      h = g - d(l)
!      d(l2:n) = d(l2:n) - h
!      f = f + h
!      !
!      !  QL transformation.
!      !
!      p = d(m)
!      c = 1.0D+00
!      c2 = c
!      el1 = e(l1)
!      s = 0.0D+00
!      mml = m - l
!      do ii = 1, mml
!         c3 = c2
!         c2 = c
!         s2 = s
!         i = m - ii
!         g = c * e(i)
!         h = c * p
!         r = pythag ( p, e(i) )
!         e(i+1) = s * r
!         s = e(i) / r
!         c = p / r
!         p = c * d(i) - s * g
!         d(i+1) = h + s * ( c * g + s * d(i) )
!         !
!         !  Form vector.
!         !
!         do k = 1, n
!            h = z(k,i+1)
!            z(k,i+1) = s * z(k,i) + c * h
!            z(k,i) = c * z(k,i) - s * h
!         end do
!      end do
!      p = - s * s2 * c3 * el1 * e(l) / dl1
!      e(l) = s * p
!      d(l) = c * p
!      tst2 = tst1 + abs ( e(l) )
!      if ( tst2 > tst1 ) then
!         go to 130
!      end if
! 220  continue
!      d(l) = d(l) + f
!   end do
!   !
!   !  Order eigenvalues and eigenvectors.
!   !
!   do ii = 2, n
!      i = ii - 1
!      k = i
!      p = d(i)
!      do j = ii, n
!         if ( d(j) < p ) then
!            k = j
!            p = d(j)
!         end if
!      end do
!      if ( k /= i ) then
!         d(k) = d(i)
!         d(i) = p
!         do j = 1, n
!            call r8_swap ( z(j,i), z(j,k) )
!         end do
!      end if
!   end do
!   return
! end subroutine tql2


! !---------------------------------------------------------------------
! ! PURPOSE: computes SQRT ( A * A + B * B ) carefully.
! !    The formula
! !    PYTHAG = sqrt ( A * A + B * B )
! !    is reasonably accurate, but can fail if, for example, A**2 is larger
! !    than the machine overflow.  The formula can lose most of its accuracy
! !    if the sum of the squares is very large or very small.
! !  Parameters:
! !    Input, real(8) :: A, B, the two legs of a right triangle.
! !    Output, real(8) :: PYTHAG, the length of the hypotenuse.
! !---------------------------------------------------------------------
! function pythag ( a, b )
!   implicit none
!   real(8) :: a
!   real(8) :: b
!   real(8) :: p
!   real(8) :: pythag
!   real(8) :: r
!   real(8) :: s
!   real(8) :: t
!   real(8) :: u
!   p = max ( abs ( a ), abs ( b ) )
!   if ( p /= 0.0D+00 ) then
!      r = ( min ( abs ( a ), abs ( b ) ) / p )**2
!      do
!         t = 4.0D+00 + r
!         if ( t == 4.0D+00 ) then
!            exit
!         end if
!         s = r / t
!         u = 1.0D+00 + 2.0D+00 * s
!         p = u * p
!         r = ( s / u )**2 * r
!      end do
!   end if
!   pythag = p
!   return
! end function pythag

! !---------------------------------------------------------------------
! ! PURPOSE: swaps two R8's.
! !  Parameters:
! !    Input/output, real(8) :: X, Y.  On output, the values of X and
! !    Y have been interchanged.
! !---------------------------------------------------------------------
! subroutine r8_swap ( x, y )
!   real(8) :: x
!   real(8) :: y
!   real(8) :: z
!   z = x
!   x = y
!   y = z
!   return
! end subroutine r8_swap



